Long-range dispersion and spatial diffusion of fault waves in the Burridge-Knopoff 

earthquake model. 
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The Burridge-Knopoff model of earthquakes has recently gained increased interest for the con- 
sistency of the predicted energy released by sismic faults, with the Gutenberg-Richter scaling law. 
The present work suggests an improvement of this model to account for long-range dispersions and 
large spatial diffusion of sismic faults. An enhancement of the threshold speed of shock waves driven 
by translated fault fronts is pointed out and shown to result from the interactions between com- 
ponents of the system situated far aways them and others. Due to the enhanced threshold speed, 
size of the sismic fault gets increased but a control effect can still be gained from tunable disper- 
sion extent irrespective of the total length of the system. To the viewpoint of the Burridge-Knopoff 
block-lattice model, this last consideration introduces the possibility of sizable but finite interactions 
among infinitely aligned massive blocks. Implications on the fault wave propagation are examined 
by numerical simulations of the improved nonlinear partial differential equation. 

PACS numbers: 91.45.Yb, 91.30.Dk, 91.30.Mv, 02.60.Lj 



Understanding natural catastrophes is one of the great- 
est challenges faced by scientists from various fields. Re- 
cently, earthquake phenomena attracted a lot of atten- 
tion to both viewpoints of the statistics of earthquake 
events Q Q and of the dynamics of sismic faults |3| . The 
first 11,1^, based on cellular automata models and involv- 
ing fractal faults, assumes a system which is in perpetual 
crisis 0, 0- The main feature of this approach is the 
prediction of large fluctuations in the avalanche sizes 4] . 
The second approach follows a deterministic description 
assuming sismic events as time evolving processes. More 
explicitely, this last approach is based on spatio-temporal 
evolutions of the sismic fault considered within equations 
which in general are of the Klein-Gordon type. The best 
known of these equations is provided by the Burridge- 
Knopoff(BK) model 6J intended for two parallel tectonic 
plates subjected to stick-slip frictions, which can slide one 
relative to the other under a velocity-weakening driving 
force. In theory, one considers the motion of massive un- 
derground blocks attached to a static surface and inter- 
acting them and others via harmonic springs of constant 
strengths. The dynamic friction force introduces nonlin- 
earity in the fault dynamics so that the corresponding 
discrete nonlinear equation writes 0, : 

Mi/„ = if({/„+i-2{/„-|-[/„_i)-fc„{/„-$({/„/T/o), (1) 

where Un is the relative sliding amplitude of the n*'' block 
of mass M , K \s the effective coupling strength felt by the 
block in the presence of two nearest-neighbour blocks, ko 
is the uniform strength of block stick to the surface and 
dots refer to time derivatives. We assume two tectonic 
plates at equilibrium, in other words we neglect resid- 
ual uniform sliding motions of blocks as they are pulled 
individually forward(or backward) relative to the upper 



plate "^l. Concentrating ourselves on the sismicity of the 
model triggered by the dynamic friction force <&([/„/ Vo), 
this force is usually written: 



(2) 
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in which Fq and Vq are two constant parameters. 
Equation exhibits rich physical properties one most 
relevant being the Gutenberg-Richter Q scaling law 
associate to the distribution of the energy released 
by earthquakes. Otherwise, previous numerical sim- 
ulations established an agreement of the resulting 
velocity profiles and wave amplitudes with real earth- 
quake processes and routes to chaos have very recently 
been investigated [Tol |. 

Discreteness is however another relevant aspect of this 
equation and constitutes the main subject of the present 
work. Discreteness connects to dispersion which is an 
intrinsic property of the system having direct incidence 
on its stability. Dispersion in the present system is as 
foundamental as it characterizes the ability of the in- 
teracting block lattice to respond to both small and 
large amplitude excitations propagating along the lat- 
tice, namely its fixes the threshold frequency and speed 
for stable excitations. This problem has been discussed in 
previous studies ^ [ij in terms of an amplification rate 
from which the group velocity was derived by linearizing 
this physical parameter about a characteristic wavector. 
As it stands, the BK model supposes each block can only 
see the two nearest-neighbour blocks which means short- 
range dispersion. Accordingly, the corresponding thresh- 
old speed will be appropriate only for short-range(SR) 
dispersive excitations. However, spring models are flex- 
ible in essence and SR descriptions often consist of ap- 
proximations of the actual physics behind these models. 
For instance, flexibility of the lattice of massive blocks 
coupled by springs can allow a given n^^ block to re- 
spond almost instantaneously to a perturbation taking 



2 



born on or crossing an m*^ block relatively far away(i.e. 
beyond nearest neighbours) from it. Numerous sismic 
events show clear evidences of an implication of such fea- 
tures and it is of crucial importance to take them into 
account. Moreover, recent advances on spring models 
make possible their theoretical descriptions in terms of 
long-range (LR) dispersions and the associate Klein- 
Gordon equation is qualitatively similar to the SR one. 
Nevertheless, the SR coupling strength K is substituted 
for a LR dispersion potential such that the system dy- 
namics now obeys the following equation: 

MUn = K{m - n){Urn - Un) - koUn - ^{Un/Vo), 

(3) 

The LR dispersion force K{m — n) is thus an adjustable 
parameter and the equation Q turns to a deformable BK 
model. One most interesting features of this deforma- 
bility is the possibility to adjust the model to various 
situations according to the strength and nature of the 
interactions among components of the system. For in- 
stance, if we take account of the block stick to the up- 
per plate and assume blocks to be relatively strongly at- 
tached to the surface, the LR dispersion function should 
quickly decrease at short distances and an exponentially 
decreasing function with increasing distance m — n can 
be an acceptable approximation. This last consideration 
is formulated assuming an expression of the form: 



r(l — r^) 



m — n. 



(4) 



where L is the spatial extent of the interaction between 
blocks. By spatial extent we understand the number of 
neighbour blocks from a given block which can be differ- 
ent from the total number of blocks forming the whole 
lattice. This precision is of particular importance here 
since we want a model that allows a control of the spa- 
tial extent of the dispersion. In addtion to L, the quan- 
tity r is another control parameter but relates instead to 
the strengths of successive couplings. For an exponential 
fall-off we need r ~ where A is a positive constant. 
Or also, we can just confine r in the interval < r < 1. 
In this way, the LR dispersion function (|3J looks quit e 
like the so-called Kac-Baker potential ifl lisl Il6|. 
Nevertheless, unlike the usual Kac-Baker potential the 
new version (0J has finite and constant magnitude for 
finite values of L. Analytically, this is traduced by the 
constraint: 



K 



(5) 



This constraint acquires particular interest in the present 
context since by fixing the total interaction between 
blocks to a finite magnitude, we avoid uncontrollable 
dispersions of the energy carried by the LR excitations 
which, otherwise, is nothing else but the total energy car- 
ried by the fault. 



In what follows we examine dispersion properties of the 
improved model (O and point out some consequences of 
the account of the LR dispersion on the fault dynamics. 
Following usual considerations •gl, 111, we linearize this 
equation assuming small-amplitude motions of blocks 
about their equilibra. These equilibra correspond to the 
most stable positions in the lattice after uniformly pulling 
the blocks, resulting in a uniform equilibrium position 
Uo = Fo/ko- Setting Un{t) = Uo + u{q) exp {naQi + fit) 
and keeping only linear terms in u and ii, equation (jS)) 
leads to the following dispersion law: 



1 ± 
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^ M ( Kl{Q) \ f2Vo 
ko\ ko [Uo 



flo = 



Fa 



Q = Qa. 



(6) 



fl in this last relation is the amplification rate for small- 
amplitude excitations and Q harmonics associate to their 
spatial dispersions. The function Kl{Q) appearing in 
is defined as: 



KLiQ) = 2j2K{i) [cos(Q^) 



(7) 



e=i 



where a is the equilibrium separation between two neigh- 
bouring blocks. Thus, a appears as a characteristic 
length scale of the model. As we will see, the function 
Kl (Q) governing spatial dispersion provides another rel- 
evant length scale. Instructively, this second length scale 
was previously considered and called stiffness length(e.g. 
denoted C in JJ ) . Below we keep the same viewpoint but 
introduce more suggestive interpretation in terms of the 
size of fault wave as it moves in a Galilean frame i.e. with 
a translated wavefront. To start let us examine the dis- 
persion law obtained in J^. The discrete sum in K{Q) 
can be calculated analytically using the identity: 



Vl{Q) = E^'^°^(^^) 



rcos(g) - _ r(-^+i)[cos((L + 1)(5) - rcos(LQ)] 



1 - 2rcos((5) + 



With help of the analytical expression of fl{Q) derived 
from this identity an explicit formulation of the disper- 
sion property of the LR system becomes trivial. In this 
goal, note the presence of two signs ± in lO which re- 
late to shock fronts moving(translating) respectively to 
the left (backward) and to the right (forward) in the block 
lattice. To clearly see these two distinct polarities in the 
dispersion properties of the fault wave, on figure Q we 
plot the two amplification rates fl± as function of the 
wavector Q/tt for the infinite-extent dispersion and for 
some arbitrary values of the LR parameter r. The zero 
dispersion mode Q = appears as lying inside a gap of 



(8) 
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FIG. 1: Dispersion laws versus reduced wavector Q/tt. From 
the upper curve in the "high-frequency" regime(upper sub- 
spectrum): r = 0.9, 0.8, 0.5, 0. 
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FIG. 2: threshold speeds as a function of r. From the lowest 
curve: L = 4, 9, 20 41, 500. 

finite widtli Afi = — il- |q-»o. The effect of this gap 
on the overall lattice dispersion is to lift the degeneracy 
of the dispersion spectrum giving rise to two separate 
sub-spectra associate to backward and forward disper- 
sion modes. However, both sub-spectra possess common 
characteristic(sound) speed which is approximately the 
slope of the linear part of their dispersion curves and 
consequently can be estimated from: 

Cl = {l/M)dQn±{Q)\Q^o, (9) 

Cl is plotted on fig. |21 as a function of r assuming dif- 
ferent values L. A stricking feature in this last figure is 
an enhancement of the threshold speed as the number 
of interacting neighbours increases. The infinite in- 
crease of the characteristic speed for infinite interaction 
extent indicates a system becoming strongly unstable as 
all components forming this system are almost sensitives 
to their mutual presence at relatively long distances them 
and others. 

Implications of this infinite increase of the characteristic 



speed become quite clear if we consider excitations with 
long wavelengths compared to the characteristic length 
scale a, then readily assumed as an ultra-violet cut-off. 
These are continuum-limit processes of the equation ^ 
and hence can be described by the continuum nonlinear 
Klein-Gordon equation: 

Un = ClU^.^u;lU-^^Un/Vo), ^o = ^- (10) 

where na x. To see the meaning of to the system 
dynamics, we perform Galilean transformation U{s, t) — > 
U{s ± vt) setting: z = 7^ = 1 ~ v^/Cl and find: 

i;. + F-i/(i±ay,) = o, ei = 4' « = 

with Y — U/Uo- Eq. reveals that £,l is precisely the 
characteristic topological size of the sismic fault, hence 
a second characteristic length scale of the model. The 
linear dependence of this quantity on Cl is suggestive 
enough as for the effect of the LR dispersion on this sec- 
ond length scale. 

To end, let us look at the consequence of an account 
of the LR dispersion on the spatial diffusion of large- 
amplitude(soliton-like) fault waves. For this purpose, the 
continuum nonlinear equation was solved (in U) ap- 
proximating the diffusion operator(second-order spatial 
derivatives) within a central-difference scheme and apply- 
ing combined fourth-order Runge-Kutta and fifth-order 
Runge-Kutta rules(so-called Runge-Kutta-Fehlberg algo- 
rithm) with error controls on time-dependent variables. 
Arbitrary values of the model parameters were chosen 
except whose dependence on the two control param- 
eters r and L were carefully treated. To this last point, 
according to fig. |21the effect of increasing the dispersion 
extent L keeping r fixed is qualitatively similar to the 
effect of an increase of r for a fixed value of L. Therefore 
results displayed below assume two different but equiv- 
alent qualitative interpretations. On figs. O 0] [S] and El 
we present some velocity profiles as the topological size 
S,L is increased by tuning either r fixing L, or L fixing r. 
Initial conditions are {V,V) = (—1.0,0). The velocity 
profile in fig. |3I can serve as a reference since it corre- 
sponds to the SR solution(r = 0, L ^ 00). The three 
next figures are spatio-temporal shapes of the velocity 
profile for r = 0.2 0.4 and 0.6 respectively (with L = 40 
nearest neighbours). As one notes, the wave carrying the 
sismic fault gradually spreads all over the space as r in- 
creases. This suggests an increasingly large diffusion of 
the fault due to a growth of the topological size of this 
wave as r increases. Therefore, numerical simulations are 
in excellent agreement with our analytical predictions. 
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FIG. 3: Velocity profile for r = 0. 



FIG. 4: Velocity profile for r = 0.2. 
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FIG. 5: Velocity profile for r = 0.4. 



FIG. 6: Velocity profile for r = 0.6. 
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